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ABSTRACT 

This paper summarizes the technical development and validation of a multiphase 
computational fluid dynamics (CFD) numerical method using the volume-of-fluid (VOF) model 
and a Lagrangian tracking model which can be employed to analyze general multiphase flow 
problems with free surface mechanism. The gas-liquid interface mass, momentum and energy 
conservations are modeled by continuum surface mechanisms. A new solution method is 
developed such that the present VOF model can be applied for all-speed flow regimes. The 
objectives of the present study are to develop and verify the fractional volume-of-fluid cell 
partitioning approach into a predictor-corrector algorithm and to demonstrate the effectiveness of 
the present innovative approach by simulating benchmark problems including the coaxial jet 
atomization. 


INTRODUCTION 

The atomization and breakup processes of liquid fuel in many space transportation and 
propulsion systems are of vital importance to combustion performance and stability. Numerical 
modeling of these flows poses a significant challenge because it requires simultaneous 
resolutions of transient liquid-gas-droplets dynamics, and the flow regimes considered can range 
from incompressible to high-speed compressible flows. The flow-process modeling is further 
complicated by surface tension, interfacial heat and mass transfer, material property variation, 
spray formation and turbulence, and their interactions. It is the aim of this study that a general 
and robust design analysis tool will be available as a result of the incorporation of advanced 
numerical and physical models in a computational fluid dynamics flow solver. 

Historically, developments of numerical methods for multiphase free surface flows were 
primarily aimed at incompressible flows utilizing the Volume of Fluid (VOF) method. [1,2]. 
These methods are often ineffective for high-speed, high-Reynolds- number flows. So far, 
multiphase free surface flows have not been considered in more general all-speed flow 
algorithms, with only a few recent exceptions [3,4,5]. In the ARICC-3D Code [3], the VOF 
method was implemented in the ALE-ICE (Arbitrary Lagran gian Eulerian-Implicit Continuous- 
fluid Eulerian method) algorithm for injector flow simulation with atomization. Due to the 
inefficiency of the ALE-ICE method, the ARICC-3D was found very time-consuming for multi- 
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element computations. Recently, the ALE-ICE method was abandoned by the code developer in 
favor of the pressure-based SIMPLE algorithm [4], In the RIPPLE computer program [5], only 
finite difference solutions to the incompressible Navier-Stokes equation are obtained on an 
Eulerian, rectangular mesh using a time-split two-step projection algorithm. Compressible fluids 
were not solved, while the surface tension is treated using a novel continuum surface force (CSF) 
model. 

Based on a preliminary study presented last year [16], a unified multi-phase numerical 
method is developed. The gas-liquid interface mass, momentum and energy conservation 
properties are modeled by continuum surface mechanisms. A new solution method is developed 
such that the present VOF model can be applied for all-speed range flows. The objectives of the 
present study are: (a) to develop and verify the fractional volume-of-fluid (VOF) cell 

partitioning approach into a predictor-corrector algorithm to deal with multiphase (gas-liquid) 
free surface flow problems; (b) to implement the developed unified algorithm in a general 
purpose computational fluid dynamics (CFD) code. Finite Difference Navier-Stokes (FDNS), 
with numerical treatment of droplet dynamics models described in [6]. 

The main feature of the present method is to combine the novel feature of the Volume of 
Fluid (VOF) method and the Eulerian/Lagrangian method used in a pressure-based particulate 
two-phase flow solver [6,9, 10] into a unified algorithm for efficient non-iterative, time- accurate 
calculations of multiphase free surface flows valid at all speeds. The proposed method 
reformulated the VOF equation to strongly couple two distinct phases (liquid and gas), and tracks 
droplets on a Lagrangian frame when spray model is required, using a unified predictor-corrector 
technique to account for the non-linear linkages through the convective contributions of VOF. 
The discontinuities within the sharp interface will be modeled as a volume force to avoid 
stiffness. Formations of droplets and tracking of droplet dynamics are handled through the same 
unified predictor-corrector procedure. Thus the new algorithm is non-iterative and is flexible in 
the handlin g of any general geometries with arbitrarily complex topology in free surfaces. The 
present method can also be applied for efficient steady state situations and can be implemented 
into any pressure-based methodology. For the purpose of algorithm development and validation, 
the proposed method will be implemented into a general purpose CFD code. Finite Difference 
Navier-Stokes FDNS [9,10], which has three dimensional and multi-zone capabilities. 


METHODOLOGY 


Governing Equation 

Conventionally, there have been some successful numerical algorithms for performing 
multiphase flow calculations involving free interfaces. In general, two basic approaches can be 
identified. One is based on the Lagrangian method and the other on the Eulerian frame work. 
The Lagrangian approach is to conform a mesh system with the sharp interface between gas and 
liquid and to track the movement of the interface according to sets of simplified equations 
governing the liquid and gas motions separately [8]. The entire Lagrangian grid network moves 
with the free surface which makes the grid generation or remeshing procedure very complex and 
CPU intensive. In the Eulerian approach, a fixed mesh system is used to resolve the sharp 
interface by the concept of a fractional volume of fluid (VOF) [11], The second approach is 
more flexible and efficient in that it allows the resolution of the free surface to become part of 
the solution. The physical processes across the sharp interface such as surface tension forces, 
droplet or wave breakup, interphase mass transfer and condensation/vaporization, etc. can be 
modeled through the concept of continuum surface force (CSF) method [7] or by the method of 
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surface reconstruction as used in the ARICC3D code [4], SOLA-VOF code [1]. The CSF 
method is adopted in the present model for its generality and computational efficiency, especially 
for 3-dimensional computations. 


Traditionally, VOF methods are mainly developed and used for low-speed flows such that 
incompressibility can be assumed. The incompressible flow assumption has limited their 
capability. To generalize, the present formulation is based on compressible flow governing 
equations. The forms of the equations are then continuously reduced to their incompressible 
forms according to the local flow conditions and the VOF solutions. This is the uniqueness of 
the present method. To illustrate this, the density-weighted averaged conservation equation of 
mass, Navier-Stokes, and scalar variables in an Eulerian frame work can be written as: 

9p m 


dt 


dXj 
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The VOF transport equation is given below. 
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where a = 1 stands for liquid and a = 0 is for gas. The interface is located at 1 > a > 0. S a 


represents the volume transfer rate across the two-phase boundaries. And, 

P m =(l-oOp,+ap, 

where p m denotes the time-mean density of the mixture, p g and pi denote gas and liquid density 
respectively, u, and u\ are the i-component of the density-weighted mean and fluemating part of 
the instantaneous velocity, <p and <J>' are the density-weighted mean and fluemating part of the 
instantan eous scalar quantities including the species concentrations, turbulence quantities and the 
gas mixture enthalpy, p is the mean pressure, S represents sources terms due to mass transfer, 
momentum transfer, species production, etc. u g represents the grid speed components used to 
simulate moving domain effects. Detailed expressions of these source terms can be foun d in 
Refs. 9 and 10. In the present study, the turbulent correlation terms, u\ u'j and u\ <()', are 


modeled by the two-equation turbulence closure models [12]. For a given solution of a field, 
equation (1) can be rewritten as the following form to maintain accurate transient from the gas 
phase flowfield to the liquid phase flow solution through the interface. 

Pin (h dp m (M - M. ) <t> 

- 1 ! — = S. , a < 0. 0 1 for compressible gas 


dt 

d(f) 

>m dt 
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Interface Resolution Model 

The numerical accuracy of the VOF method depends highly on the interface resolution. To 
prevent the solution from becoming too smearing due to numerical diffusion, a compression 
procedure is developed to perform VOF interface rescaling such that the average thickness 
within the interface (0.1 < a < 0.9) is kept constant through out the computation. The interface a 
solution compression procedure is expressed as: 

= Afax{o, Mm[l, 0.5+/ (a oW - 0.5)]} 
and 

{Interface Thickness) new 
{Interface Thickness) ^ 


General Continuum Surface Tension Force Model 

The surface tension forces in the continuum surface force model are formulated as 
continuous body forces across the interface. These forces can be written as, 


F, = -a Vn 


F„ = -o V n 


A 


a. 


a y + 


a. 


V ' Jfo r 2D axisymetric only 


F ’ = -a Vn 


a z , for 3D case only 


where a = surface tension coefficient and V n = a« + a w + a 


Interface Mass Transfer Model 

To simulate the mass transfer effects along the interface of the gas and liquid phases, the 
source term of Eq. (2) is modeled. The sign and magnitude of the VOF source term, S a in Eq. 
(2), depend on the physical processes involved. For spray atomization applications, for instance, 
it would represent the volume stripping rate which can be modeled by atomization correlations. 
For spray coating, condensation and chemical vapor deposition processes, on the other hand, 
positive volume flow rates would result. 

There are basically two ways of modeling the VOF source term. One method is to treat the 
inter-phase volume flow rate as convection process along the interface. That is, 

s a = y; • Voc 

where V s denotes the volume exchange velocity. This method allows the source term and the 
convection terms be combined so that the source effects are handled implicitly. However, the 
accuracy of this approach depends entirely on the assignment of the volume exchange velocity. 
Vectors directly normal to the interface do not always guarantee good solution. 

The second approach, which is the most straight forward treatment, solve the VOF equation 
with the source term specified explicitly. This approach is general and provides good numerical 
accuracy. The only drawback of this method is the possible numerical instability related to high 
volume flow rates. Time step size must be small enough in order to solve the VOF equation with 
large source term. 
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Numerical Algorithm 
Time-Marching Scheme 


For simplicity, the density-weighted conservation equation of mass, momentum, enthalpy, 
and scalar variables in an Eulerian coordinate can be written as follows: 



dp<t> 

dt 


JT {pU,) = s ’ 

+ ^-(p £/,♦) = 5 


♦ 


where consists of the diffusion fluxes and the source terms. Using general 0 method for 
temporal discretization, the above transport equations are discretized as: 


( P "*‘-p")/Ax+e ( P f/)”‘+(i-eXp(/)>e5”'+(i-e)5; 

[w * 1 -<p*)"]/ A»+e (pt/c 1 +(i-e)(p£/ 4 >): =» s :" +(1-9)5; 


where 0 = 1/2 represent the time-centered Crank-Nicholson scheme which is second-order 
accurate in time. This scheme is used in the present study. By substituting the continuity 
equation into the scalar transport equation, the delta form of the advection equation of the scalar 
transport equation is obtained: 


f-(i-8)( P t/); 

_Af 

= “[(pH )" + A(p(/ +(l-8)A(pC/)<|>, +8 Sf*' 
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The predictor and corrector of the above equation are then written as: 
i£-(l-e)(p£/); + 0 [(pI/)" + A(pt/)]A«, 

= -[(p u) m + A (pc/ )]<t» x + (1-0) A(p I/)*, + 5 ; 


and 


H^-(i-e)(pc/)“ 


4>’+e [(pur+Alpc/^^es’, 


In the corrector stage, pressure field, velocity vectors, density field and enthalpy are updated by 
solving the pressure correction equation and the scalar correction equations. For wave 
propagation problems, 3 to 4 correctors are sufficient for convergence. 


High- Order TVD Scheme 

In the present study all test cases were analyze with a high-order upwind TVD Scheme. 
Only the convection terms are modeled using the TVD flux limiters. The convection terms of the 
governing equations can be expressed by finite difference approximation as: 

At = fi+v 2 ~ fi-m fy+v 2 — ^1-1/2 
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where f and h represent first-order fluxes and TVD flux limiters respectively. The TVD flux 
limiters are functioned as anti-diffusion terms to recover the scheme to high-order accuracy. The 
first-order fluxes and the TVD flux limiters are given below. 


f 1 + 1/2 = max {0,(pC/ }<!>,. + max 


h 


1 + 1/2 




{o.-ipV),.,, 

+ cx (^<|> l+ i/2 — d§i-u2 )} > U — 0 
+ CL (^<)), + 1 / 2 — </<!>, + 3/2 )}* U <0 


where the minmod functions in the TVD flux limiters are written as: 


dtf+m = sign (A0, +1/2 )niax{o,min[|A(|), +1/2 |,p sign (A<|), +1/2 )A0, +1/2t1 ]} 


The order of accuracy of this scheme is determined by the parameters a and p. Only the second- 
order and third-order upwind schemes were used in this study. That is, 

f-1, 2 nd- order upwind 


a = < 


1 


+— , 3rd -order upwind 

3— a 

6 = Compression — Factor— 

1-a 


The compression factor, p, is used to sharpen the contact discontinuities and slip streams for 
better wave tracking. 


VALIDATIONS 


Microgravity Fluid Tank Test Cases 

To validate the present VOF method within the compressible and incompressible flow 
solver, some fluid tank benchmark test cases were investigated and compared with close form 
solutions. The key mechanisms involved here are the surface tension forces and the effects of 
wall contact angle. The liquid surface tension forces are modeled using the proposed continuum 
surface force (CSF) model and the wall contact angles are provided by the material properties 
and through the use of Young’s equation. 

Rotating Ethanol Tank : 

In this test case, a rotating cylindrical ethanol tank (46.95% full) with radius of 3 cm and 
height of 2 cm was simulated. The rotational speed of the tank is 10 rpm. The surface tension 
coefficient of ethanol and the level of gravity force (g/go) are 2.28xl0' 2 and 10‘ 3 N/m 
respectively. A grid size of 31x31 with near wall grid clustering was employed for the 
computation. 2,000 time steps with a time step size of 2.0 were required for the surface waves to 
subside and reach a steady state solution. Figure 1 shows the predicted flow field and liquid 
surface shape and comparison with close form solution. Good agreement is reveal in the data 
comparison. 
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Figure 1. Rotating ethanol tank: 10' 3 g 0 andlOrpm. 

Rotating Liquid Helium Tank : 

In this test case, a rotating cylindrical helium tank (50% full) with inner and outer radii of 
18.48 cm and 77.48 cm respectively, and height of 136 cm was investigated. The rotational 
speed of the tank is 1 rpm. The surface tension coefficient of helium (normal fluid) at 1.5 degree 
K and the level of gravity force (g/go) are 3.319xl0‘ 4 N/m and 10" 6 respectively. A grid size of 
31x31 with near wall grid clustering was employed for the computation. 3,000 time steps with a 
time step size of 5.0 were required for the surface waves to subside and reach a steady state 
solution. Figure 2 shows the predicted flow field and liquid surface shape and comparison with 
close form solution. Good agreement is reveal in the data comparison. 




Figure 2. Rotating L-He tank: 10' 6 g 0 and 1 rpm. 
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Cryogenic Liquid Helium Fluid Tank : 

The working fluid and fluid properties of this test case is identical to that of the case 2 above. 
A dewar helium tank (with inner and outer radii of 12 cm and 70 cm respectively and height of 
150 cm) under 10" 7 go and 0.1 rpm operating condition was tested. This case also provides a 
validation test of the surface tension force model for non-orthogonal grid systems. A steady state 
solution was obtained in 3,000 time steps with a time step size of 5.0. Figure 3 shows the 
predicted liquid surface shape and flow field. This solution is in good agreement with the 
numerical solution presented in Ref. 14. 




Figure 3. Rotating dewar solution: 10' 7 g 0 and 0.1 rpm. 

Rising Air Bubble in a Tank : 

In this simulation, a spherical air bubble with radius of 20 cm was released at the center of a 
tank under lg 0 room temperature conditions. A 31x31 grid with time step size of 0. 1 was used in 
the computation. This case is used to check the correctness of the pressure solutions across the 
interface and inside the bubble. Results of the predicted liquid interface, flowfield and pressure 
field for time level of 0.5 sec and 0.7 sec are given in Figure 4. It is clear that the predicted 
pressure inside the bubble stay almost constant and the shape of the pressure contour follows that 
of the bubble boundary. This simulation is considered to be physically correct. 
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Figure 4. Rising air bubble simulation: lg 0 conditions. 




Coaxial Liquid Jet Atomization 

A cold flow case of coaxial liquid jet atomization reported by Liang and Schuman [15] is 
employed to test the current numerical methodology. Figure 5 shows the injector configuration. 
We use our unstructured grid code to simulate this multi-flow passage problem and expect to 
simulate multi-injector combustion applications by taking its advantage of grids flexibility for 
complex combustor geometry. In the current stage, the breakup rate is fixed as 10% of the liquid 
jet velocity and the droplet Sauter mean diameter is 100pm with Rosin- Rammler distribution. 
This restriction can be removed by using the wave instability model of Reitz et al. [17]. 

H 2 0 and C0 2 are used as the liquid and gas in this study with injection velocity 12.14 
m/s and 166.4 m/s respectively. 100 time steps are executed without injecting of particle but 
with liquid jet breakup. After that we ran another 50 time steps with particle injection from the 
surface and smaller time-step size (20ps). Figure 6 shows the velocity vectors and liquid core. 
Liquid jet velocity decreases due to the nozzle expansion and then is accelerated by the large slip 
velocity of gas phase. The numerical model is stable even for large time-step size and 
incompressible (liquid) /compressible (gas) flow situation. Figure 7 presents the particle 
trajectories at time = 1 ms. Particle turbulent dispersion effect has been included and large 
particle dispersion is shown due to high turbulent intensity of gas phase. The numerical 
prediction looks reasonable for the main flow feature. Further comparisons with benchmark data 
from Penn. State University are currently underway. 
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Figure 7. Particle trajectories at t = 1 ms. 


10 





CONCLUSIONS 


This paper summarizes the technical development and validation of a unified multiphase 
computational fluid dynamics (CFD) numerical method using volume-of-fluid (VOF) model 
which can be employed to analyze general multiphase flow problems with free surface 
mechanism. A time-accurate multiphase viscous flow solution method has been developed and 
validated with a wide range of benchmark test cases including the coaxial jet atomization. The 
gas-liquid interface mass, momentum and energy conservation properties are modeled by 
continuum surface mechanisms. A new solution method is developed such that the present VOF 
model can be applied for all-speed flow regimes. The objectives of the present study are to 
develop and verify the fractional volume-of-fluid cell partitioning approach into a predictor- 
corrector algorithm to deal with multiphase (gas-liquid) free surface flow problems; and to 
couple this unified algorithm with a non-iterative Lagrangian model. The successful 
implementation of this study will provide a reliable numerical model for the multiphase flow 
analysis as well as the fundamental building block for a analytical tool which may lead us toward 
the understanding of complex multiphase flow physics that is sometimes hard to visualize 
experimentally. 
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